#pragma once
#include "cem/constants.h"
#include <array>
#include <libmesh/vector_value.h>

namespace cem {
namespace material {

// 工具函数: 角度转弧度
inline double deg_to_rad(double deg) {
    return deg * cem::constants::pi / 180.0;
}

// 旋转矩阵 R = Rz(strike) * Rx(dip) * Ry(slant)
std::array<std::array<double, 3>, 3>
build_rotation_matrix(double strike_deg, double dip_deg, double slant_deg);

// T = R * diag(p1,p2,p3) * R^T
std::array<std::array<double, 3>, 3>
build_tensor_from_principal(double p1, double p2, double p3,
                            double strike_deg, double dip_deg, double slant_deg);

// T^{-1} = R * diag(1/p1,1/p2,1/p3) * R^T（带安全倒数）
std::array<std::array<double, 3>, 3>
build_inv_tensor_from_principal(double p1, double p2, double p3,
                                double strike_deg, double dip_deg, double slant_deg);

// 收缩: g_i^T M g_j
double contract_gradients(const libMesh::RealGradient & gi,
                          const libMesh::RealGradient & gj,
                          const std::array<std::array<double, 3>, 3> & M,
                          unsigned int dim);

// 收缩: v_i^T M v_j
double contract_vectors(const libMesh::RealVectorValue & vi,
                        const libMesh::RealVectorValue & vj,
                        const std::array<std::array<double, 3>, 3> & M,
                        unsigned int dim);

} // namespace material
} // namespace cem